%Structual estimation (main file),including variable cost savings from the
%dominant local vendor



clear
tic

choice=1:13;
beta=0.95;
T=100;
length=0:T;
df=beta.^length;
t0=2006;
tN=2009;
tbar=tN-t0+1;
%AI_fs=[51:63,65:66,68:72,76:94,96:119,121:145,147:150];


start=1;
page=start:(start+29);
narray=numel(page);
%narray=s;
%coef=zeros(narray,2);
%ste=zeros(narray,2);
%fv=zeros(narray,1);

fs0=dlmread('fs1.csv');
fs0=sortrows(fs0,[1,2,5]);
ddom0=dlmread('mks1.csv');
ddom0=sortrows(ddom0,[1,2,3]);
FSTA=zeros(size(fs0,1),size(fs0,2),narray);
DOM=zeros(size(ddom0,1),size(ddom0,2),narray);


% Construct other inputs from "regvar*"
bd=fs0(:,3,1);  %bdtot: 25%: 38; 50%: 89; 75%: 200
%nobd=bd(:,ones(T+1,1));
%bdtype=(bd<=38&bd>0)+(bd<=89&bd>38)*2+(bd<=200&bd>89)*3+(bd>200)*4;
bdtype=(bd<=25&bd>0)+(bd<=80&bd>25)*2+(bd<=193&bd>80)*3+(bd>193)*4; %bdtot: 25%: 25; 50%: 80; 75%: 193

%### UPDATED on 05/07/2020:Calculate the mktdom at the beginning here NOT using the one in the policy
%function because we need the endogenous mktdom (including own choice) so
%each hospital in the same market has the same mktdom
textFileName2006=sprintf('regvar6_sing.raw');  %the input file     
textFileName2007=sprintf('regvar7_sing.raw');  %the input file     
textFileName2008=sprintf('regvar8_sing.raw');  %the input file    
textFileName2009=sprintf('regvar9_sing.raw');  %the input file    


char2006=dlmread(textFileName2006);   %the matrix of characteristics of stand-alone hospitals: ahaid year bdtot prechs adj_vnd hsa 
char2007=dlmread(textFileName2007);
char2008=dlmread(textFileName2008);
char2009=dlmread(textFileName2009);

allchar=[char2006;char2007;char2008;char2009];
allchar=sortrows(allchar,[1,2]); 
pre=allchar(:,4);
act=allchar(:,5);  %extract the current action to compute the probability of the actual choice

nowDOM=zeros(size(allchar,1),numel(choice));
for i=1:size(allchar,1)
    tmpmkt=allchar(allchar(:,end)==allchar(allchar(:,1)==allchar(i,1)&allchar(:,2)==allchar(i,2),end),4);
    tmpmkt(tmpmkt==1)=[];
    tmpmkt(tmpmkt==2)=[];
    tmpmkt(tmpmkt==13)=[];
    if isempty(tmpmkt)==0
       occur=histc(tmpmkt,choice)';
       occur=reshape(occur,1,[]);           
       mktdom=choice(occur==max(occur));
       nowDOM(i,mktdom)=1;
    end    
end
nowDOM=[allchar(:,1:2),nowDOM];

% Define cuboid for fminsearch
df_cub=repmat(df,size(FSTA,1),1,narray);
bd_cub=repmat(bd,1,T+1,narray);
bdtype_cub=repmat(bdtype,1,T+1,narray);
gg_cub=zeros(size(DOM,1),T+1,narray);
fchs_cub=zeros(size(FSTA,1),T+1,narray);
ind_cub=zeros(size(FSTA,1),T+1,narray);
idx_sw_cub=zeros(size(FSTA,1),T+1,narray);
eerr_cub=zeros(size(FSTA,1),T+1,narray);
curdom_cub=zeros(size(FSTA,1),T+1,narray); %### UPDATED on 11/01/2020: for variable cost savings


%### UPDATED on 07/14/2020: create variables for market (un)observables
%textFileName_mktvar=sprintf('mktvar.raw');  %the input file 
textFileName_mktvar=sprintf('mktvar_pophhi.raw');  %the input file
mktvar=dlmread(textFileName_mktvar); %the matrix: ahaid year hhi mktcat totpop65
mktcat_tmp=zeros(size(fs0,1),1);
hhi_tmp=zeros(size(fs0,1),1);
pop65_tmp=zeros(size(fs0,1),1);
for i=1:size(mktvar,1)
    hhi_tmp(fs0(:,1)==mktvar(i,1)&fs0(:,2)==mktvar(i,2),:)=mktvar(i,3);
    mktcat_tmp(fs0(:,1)==mktvar(i,1)&fs0(:,2)==mktvar(i,2),:)=mktvar(i,4);
    pop65_tmp(fs0(:,1)==mktvar(i,1)&fs0(:,2)==mktvar(i,2),:)=mktvar(i,5);
end
hhi_tmp(fs0(:,3)==0)=0;
mktcat_tmp(fs0(:,3)==0)=0;
pop65_tmp(fs0(:,3)==0)=0;
%hhi=hhi_tmp(:,ones(101,1));
%mktcat=mktcat_tmp(:,ones(101,1));
hhi_cub=repmat(hhi_tmp,1,T+1,narray);
mktcat_cub=repmat(mktcat_tmp,1,T+1,narray);
pop65_cub=repmat(pop65_tmp,1,T+1,narray);


% ### UPDATED on 11/01/2020: incorporate the concurrent mktdom to calculate
% the variable cost savings for mktdom vendors
CURDOM=zeros(size(ddom0,1),size(ddom0,2),narray);


%### UPDATED on 05/07/2020: construct the column with the mktdom of the
%current hospital
A=[kron(nowDOM(:,1:2),ones(numel(choice),1)),repmat(choice',size(nowDOM,1),1),reshape(nowDOM(:,3:end)',[],1)];  %matrix A: [ahaid year vid mktdom]
%{
actprob=zeros(size(eer0,1),1);
for k=1:size(eer0,1)
    if eer0(k,3)~=0
       actprob(k)=A(A(:,1)==eer0(k,1)&A(:,2)==eer0(k,2)&A(:,3)==eer0(k,3),end);
    end
end
MKTDOM=[repmat(actprob,1,tbar-1),zeros(size(ddom0,1),T-tbar+2)];
%}
MKTDOM=sortrows(A,[1,2,3]);
MKTDOM(ddom0(:,3)==0,3)=0;
MKTDOM_cub=repmat(MKTDOM(:,end),1,T+1,narray);
clear A


eer0=dlmread('eerr1.csv'); %matrix: ahaid year option_given pr1 pr2 ... pr100
EERR=zeros(size(eer0,1),size(eer0,2)+1,narray);   %#### UPDATED on 05/05/2017: (Take log!!!) the output for the expected logit error: ahaid, option given, prob@t=1, prob@t=2, ..., prob@t=100
clear fs0 ddom0 eer0 
clear hhi_tmp mktcat_tmp pop65_tmp
for c=1:narray
    s=page(c);
    textFileName1=sprintf('fs%d.csv',s);
    textFileName2=sprintf('mks%d.csv',s);   
    textFileName3=sprintf('eerr%d.csv',s);
    textFileName_curdom=sprintf('ddom%d.csv',s);  % ### UPDATED on 11/01/2020: incorporate the concurrent mktdom
   
    tempEERR=dlmread(textFileName3);
    tempEERR(tempEERR==0)=1;
    EERR(:,:,c)=[tempEERR(:,1:3),zeros(size(tempEERR,1),1),log(tempEERR(:,4:end))]; %#### UPDATED on 03/13/2017:include one zero column before all the non-zero expected errors because the first one flow has no expected error.
    EERR(:,:,c)=sortrows(EERR(:,:,c),[1,2,3]);    
    FSTA(:,:,c)=dlmread(textFileName1);  %the matrix: ahaid year bdtot prechs action_given fs1 fs2...fs100
    DOM(:,:,c)=dlmread(textFileName2);  %the matrix: ahaid year action_given ddom_pre ddom0 ddom1 ddom2...ddom99 
    CURDOM(:,:,c)=dlmread(textFileName_curdom);  % ### UPDATED on 11/01/2020: %the matrix: ahaid year action_given ddom_pre ddom0 ddom1 ddom2...ddom99 
    FSTA(:,:,c)=sortrows(FSTA(:,:,c),[1,2,5]);
    DOM(:,:,c)=sortrows(DOM(:,:,c),[1,2,3]);
    CURDOM(:,:,c)=sortrows(CURDOM(:,:,c),[1,2,3]);
    
    %### UPDATED on 06/22/2020: try contructing cubic so as to remove the for loops in fminsearch 
    % ev of mktdom
    gg_cub(:,:,c)=DOM(:,4:end,c);
    % cost-saving per bed
    fchs_cub(:,:,c)=FSTA(:,4:end-1,c)+11;
    % sunk costs
    fir=FSTA(:,4:end-1,c);
    lat=FSTA(:,5:end,c);
    ind_cub(:,:,c)=(fir~=lat).*lat;
    % switching cost
    idx_sw_cub(:,:,c)=(fir~=lat).*(fir>1);
    % the expected errors
    eerr_cub(:,:,c)=EERR(:,4:end,c);  
    % ### UDPATED on 11/01/2020: variable cost savings
    curdom_cub(:,:,c)=CURDOM(:,4:end,c);
end
clear fir lat

%hbed=FSTA(:,3,1); %bdtot: 25%: 25; 50%: 80; 75%: 193
%bdtype=(hbed<=25&hbed>0)+(hbed<=80&hbed>25)*2+(hbed<=193&hbed>80)*3+(hbed>193)*4;


%### UPDATED on 05/07/2020: move this part to the fron to calculate
%absolute mktdom
%{
textFileName26=sprintf('regvar6_sing.raw');  %the input file  
textFileName27=sprintf('regvar7_sing.raw');  %the input file     
textFileName28=sprintf('regvar8_sing.raw');  %the input file     
textFileName29=sprintf('regvar9_sing.raw');  %the input file    

char6=dlmread(textFileName26); %the matrix of characteristics of stand-alone hospitals: ahaid year bdtot prechs adj_vnd hsa
char7=dlmread(textFileName27);    
char8=dlmread(textFileName28);
char9=dlmread(textFileName29);
%}


%%%%%%%%%%%%%%%%%%%%%%%%%%
%with 26 parameters
%theta0=[-6.752912 -6.311207 -3.715860 -4.035594 -6.650967 -6.566256 -7.581982 -5.019611 -5.121039 -5.779016 -4.164671 -4.502513 0.000369 0.000534 0.000058 -0.000058 0.000528 0.000580 0.000505 0.000105 0.000477 0.000506 0.000357 0.000294 0.062659 -1.8];
%theta0=[-5.869513 -5.330296 -2.948035 -3.377029 -5.751234 -5.651111 -6.713429 -4.106970 -4.216721 -4.870577 -3.188701 -3.664923 0.000197 0.000399 -0.000042 -0.000093 0.000394 0.000460 0.000374 -0.000148 0.000345 0.000378 0.000253 0.000127 0.028334 -1.428617];
 

 %theta0=[-6.0418682	-5.7275919	-4.2186812	-3.8558114	-5.7440653	-5.7927766 ...
 % -6.3105595	-4.9054132	-4.9710925	-5.2413514	-4.2134675	-3.7461579 ...
 % 0.000326417	0.000379318	-7.42E-05	-0.001186604	0.000380924	0.000458992 ...
 % 0.000466934	-0.000582455	0.000338908	0.000397496	0.00019452	-5.66E-05 ... 
 % 0.0131	  -0.018  0.027  	0.00043 	-0.0055  	-0.0073 ...	% UPDATED on 05/05/2020: add VARIOUS decreasing sunk cost (2-7)
 % -0.0093	0.0041 	0.00061	 -0.0053	0.0051	 0.000178 ...   % UPDATED on 05/05/2020: add VARIOUS decreasing sunk cost (8-13)
 % -0.0013 ...   %### UPDATED on 05/06/2020: add interaction of desunk and mktdom
 % 0.21580906	0.12910263	0.083096889	0.013831496	-1.0373199	-1.8934718	-0.98035968	-0.27365732];  %### UPDATED on 05/05/2020: add decreasing sunk cost

 
 theta0=[-5.3377957	-5.4169709	-5.0456836	-4.5787161	-4.9396757	-5.7800033 ...
  -5.1499117	-5.4216185	-3.7909583	-4.0971877	-4.5513865	-2.5642657 ...
  0.000723144	0.001076152	-0.000227134	-0.000806689	0.001165173	0.001066544 ...
  0.001304364	0.000100488	0.000950094	0.000983083	0.000685736	0.000726507 ... 
  -0.070665336	-0.015375363	0.13642023	0.13677074	-0.036445829	0.029318287 ...	% ### UPDATED on 07/14/2020: coefficients for HHI
-0.17451207	0.098704371	0.048806374	-0.1073208	0.090956289	-0.037909795 ...	
-0.003516718	-0.002002271	0.012521031	0.008135619	-0.010756333	-0.00106592 ...	% ### UPDATED on 07/14/2020: coefficients for totpop65
-0.017805938	0.00778993	-0.006986701	-0.005850727	0.004583828	-0.014687173 ...	
0.18299088	0.5215033	2.5220589 ...             % ### UPDATED on 07/14/2020: coefficients for market categories
   	-0.000137714	-0.000596333	0.000626649	3.56E-05	0.000131689 ... % ### UPDATED on 12/19/2020: mktdomXvendorFEs (2-7); for vendor 2 and 13 (Self-developed and Others): should be ZERO)
   -0.000497888	-6.05E-05	0.000126601	0.000244057	7.71E-05	 ...  % ### UPDATED on 12/19/2020: mktdomXvendorFEs (8-13)
  0.066132757	0.076975354	0.10207997	0.06776215	-1.4009128	-1.7142822	-1.485853	-1.2361045];  %### UPDATED on 05/05/2020: add decreasing sunk cost



 %theta0=[-5.6727427	-5.4148654	-2.8961247	-3.2231362	-5.8298762	-5.7155734 ...
 % -6.9408619	-4.1749532	-4.2832556	-4.8858257	-3.2607835	-3.4358166 ...
 % 0.00195965	0.00363447	-0.00157124	-0.00305057	0.00370564	0.00420194 ...
 % 0.00413653	-0.001424	0.00309183	0.00308829	0.00197276	0.00138111 ... 
 % 0.065  0.055  0.045  0.035  -1.6  -1.5  -1.4  -1.3];  %### UPDATED on 08/11/2018: the last eight: mktdom (4-1) and swc (1-4)


%theta0=[-4.9:-0.1:-6,0.0008:0.0001:0.0019,0.1];

%options = optimset('MaxIter',1e+10000000000000000000000,'MaxFunEvals',1e+10000000000000000000);
options = optimoptions(@fminunc,'MaxIter',1e+10000000000000000000000,'MaxFunEvals',1e+10000000000000000000,'Algorithm','quasi-newton');
%[theta,fval] = fminsearch(@(theta) logL_ReSTAT_vcdom_fullmkt(gg_cub,bd_cub,bdtype_cub,fchs_cub,ind_cub,idx_sw_cub,MKTDOM_cub,mktcat_cub,hhi_cub,pop65_cub,eerr_cub,df_cub,pre,act,theta),theta0,options);
[theta,fval] = fminunc(@(theta) logL_ReSTAT_VARvcsav(gg_cub,bd_cub,bdtype_cub,fchs_cub,ind_cub,idx_sw_cub,curdom_cub,mktcat_cub,hhi_cub,pop65_cub,eerr_cub,df_cub,pre,act,theta),theta0,options);


textFileName_theta=sprintf('theta_VARvcsavFmkt_new_w5fminunc%d.csv',start);  
dlmwrite(textFileName_theta,[theta fval],'delimiter', ',', 'precision', 8)


%%%%%%%%%%%%%%%%%%
%Compute the standard error
epsilon=1e-6;  %define the disturbance
pr=prob_ReSTAT_VARvcsav(gg_cub,bd_cub,bdtype_cub,fchs_cub,ind_cub,idx_sw_cub,curdom_cub,mktcat_cub,hhi_cub,pop65_cub,eerr_cub,df_cub,pre,act,theta); %pr26 is the probability function for 26 parameters.
oldpr=pr;
oldlnp=log(pr); %the probability given the original estimates
oldq=theta; %the original estimates
dldq=[];
for i=1:size(theta,2)   %the following for-loop computes the 1st derivative w.r.t. theta
    theta=oldq;
    theta(i)=oldq(i)+epsilon;
    pr1=prob_ReSTAT_VARvcsav(gg_cub,bd_cub,bdtype_cub,fchs_cub,ind_cub,idx_sw_cub,curdom_cub,mktcat_cub,hhi_cub,pop65_cub,eerr_cub,df_cub,pre,act,theta);
    newlnp=log(pr1);
    dldq=[dldq (newlnp-oldlnp)/epsilon];
end
%A=dldq'*dldq;
%var_q=diag(inv(dldq'*dldq),0);
se_q=sqrt(diag(inv(dldq'*dldq),0));   %the standard error is the sum of the outer product of first derivative

textFileName_se=sprintf('se_VARvcsavFmkt_new_w5fminunc%d.csv',start);  %the input file  
dlmwrite(textFileName_se,se_q,'delimiter', ',', 'precision', 8)
%}

toc
%diary off

